home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgelss.f < prev    next >
Text File  |  1996-07-19  |  21KB  |  604 lines
  1.       SUBROUTINE DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
  2.      $                   WORK, LWORK, INFO )
  3. *
  4. *  -- LAPACK driver routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
  11.       DOUBLE PRECISION   RCOND
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DGELSS computes the minimum norm solution to a real linear least
  21. *  squares problem:
  22. *
  23. *  Minimize 2-norm(| b - A*x |).
  24. *
  25. *  using the singular value decomposition (SVD) of A. A is an M-by-N
  26. *  matrix which may be rank-deficient.
  27. *
  28. *  Several right hand side vectors b and solution vectors x can be
  29. *  handled in a single call; they are stored as the columns of the
  30. *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
  31. *  X.
  32. *
  33. *  The effective rank of A is determined by treating as zero those
  34. *  singular values which are less than RCOND times the largest singular
  35. *  value.
  36. *
  37. *  Arguments
  38. *  =========
  39. *
  40. *  M       (input) INTEGER
  41. *          The number of rows of the matrix A. M >= 0.
  42. *
  43. *  N       (input) INTEGER
  44. *          The number of columns of the matrix A. N >= 0.
  45. *
  46. *  NRHS    (input) INTEGER
  47. *          The number of right hand sides, i.e., the number of columns
  48. *          of the matrices B and X. NRHS >= 0.
  49. *
  50. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  51. *          On entry, the M-by-N matrix A.
  52. *          On exit, the first min(m,n) rows of A are overwritten with
  53. *          its right singular vectors, stored rowwise.
  54. *
  55. *  LDA     (input) INTEGER
  56. *          The leading dimension of the array A.  LDA >= max(1,M).
  57. *
  58. *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
  59. *          On entry, the M-by-NRHS right hand side matrix B.
  60. *          On exit, B is overwritten by the N-by-NRHS solution
  61. *          matrix X.  If m >= n and RANK = n, the residual
  62. *          sum-of-squares for the solution in the i-th column is given
  63. *          by the sum of squares of elements n+1:m in that column.
  64. *
  65. *  LDB     (input) INTEGER
  66. *          The leading dimension of the array B. LDB >= max(1,max(M,N)).
  67. *
  68. *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
  69. *          The singular values of A in decreasing order.
  70. *          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
  71. *
  72. *  RCOND   (input) DOUBLE PRECISION
  73. *          RCOND is used to determine the effective rank of A.
  74. *          Singular values S(i) <= RCOND*S(1) are treated as zero.
  75. *          If RCOND < 0, machine precision is used instead.
  76. *
  77. *  RANK    (output) INTEGER
  78. *          The effective rank of A, i.e., the number of singular values
  79. *          which are greater than RCOND*S(1).
  80. *
  81. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  82. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  83. *
  84. *  LWORK   (input) INTEGER
  85. *          The dimension of the array WORK. LWORK >= 1, and also:
  86. *          LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
  87. *          For good performance, LWORK should generally be larger.
  88. *
  89. *  INFO    (output) INTEGER
  90. *          = 0:  successful exit
  91. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  92. *          > 0:  the algorithm for computing the SVD failed to converge;
  93. *                if INFO = i, i off-diagonal elements of an intermediate
  94. *                bidiagonal form did not converge to zero.
  95. *
  96. *  =====================================================================
  97. *
  98. *     .. Parameters ..
  99.       DOUBLE PRECISION   ZERO, ONE
  100.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  101. *     ..
  102. *     .. Local Scalars ..
  103.       INTEGER            BDSPAC, BL, CHUNK, I, IASCL, IBSCL, IE, IL,
  104.      $                   ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
  105.      $                   MAXWRK, MINMN, MINWRK, MM, MNTHR
  106.       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
  107. *     ..
  108. *     .. Local Arrays ..
  109.       DOUBLE PRECISION   VDUM( 1 )
  110. *     ..
  111. *     .. External Subroutines ..
  112.       EXTERNAL           DBDSQR, DCOPY, DGEBRD, DGELQF, DGEMM, DGEMV,
  113.      $                   DGEQRF, DLABAD, DLACPY, DLASCL, DLASET, DORGBR,
  114.      $                   DORMBR, DORMLQ, DORMQR, DRSCL, XERBLA
  115. *     ..
  116. *     .. External Functions ..
  117.       INTEGER            ILAENV
  118.       DOUBLE PRECISION   DLAMCH, DLANGE
  119.       EXTERNAL           ILAENV, DLAMCH, DLANGE
  120. *     ..
  121. *     .. Intrinsic Functions ..
  122.       INTRINSIC          MAX, MIN
  123. *     ..
  124. *     .. Executable Statements ..
  125. *
  126. *     Test the input arguments
  127. *
  128.       INFO = 0
  129.       MINMN = MIN( M, N )
  130.       MAXMN = MAX( M, N )
  131.       MNTHR = ILAENV( 6, 'DGELSS', ' ', M, N, NRHS, -1 )
  132.       IF( M.LT.0 ) THEN
  133.          INFO = -1
  134.       ELSE IF( N.LT.0 ) THEN
  135.          INFO = -2
  136.       ELSE IF( NRHS.LT.0 ) THEN
  137.          INFO = -3
  138.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  139.          INFO = -5
  140.       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
  141.          INFO = -7
  142.       END IF
  143. *
  144. *     Compute workspace
  145. *      (Note: Comments in the code beginning "Workspace:" describe the
  146. *       minimal amount of workspace needed at that point in the code,
  147. *       as well as the preferred amount for good performance.
  148. *       NB refers to the optimal block size for the immediately
  149. *       following subroutine, as returned by ILAENV.)
  150. *
  151.       MINWRK = 1
  152.       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
  153.          MAXWRK = 0
  154.          MM = M
  155.          IF( M.GE.N .AND. M.GE.MNTHR ) THEN
  156. *
  157. *           Path 1a - overdetermined, with many more rows than columns
  158. *
  159.             MM = N
  160.             MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N,
  161.      $               -1, -1 ) )
  162.             MAXWRK = MAX( MAXWRK, N+NRHS*
  163.      $               ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) )
  164.          END IF
  165.          IF( M.GE.N ) THEN
  166. *
  167. *           Path 1 - overdetermined or exactly determined
  168. *
  169. *           Compute workspace neede for DBDSQR
  170. *
  171.             BDSPAC = MAX( 1, 5*N-4 )
  172.             MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
  173.      $               ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) )
  174.             MAXWRK = MAX( MAXWRK, 3*N+NRHS*
  175.      $               ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) )
  176.             MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
  177.      $               ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  178.             MAXWRK = MAX( MAXWRK, BDSPAC )
  179.             MAXWRK = MAX( MAXWRK, N*NRHS )
  180.             MINWRK = MAX( 3*N+MM, 3*N+NRHS, BDSPAC )
  181.             MAXWRK = MAX( MINWRK, MAXWRK )
  182.          END IF
  183.          IF( N.GT.M ) THEN
  184. *
  185. *           Compute workspace neede for DBDSQR
  186. *
  187.             BDSPAC = MAX( 1, 5*M-4 )
  188.             MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC )
  189.             IF( N.GE.MNTHR ) THEN
  190. *
  191. *              Path 2a - underdetermined, with many more columns
  192. *              than rows
  193. *
  194.                MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  195.                MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
  196.      $                  ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  197.                MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
  198.      $                  ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) )
  199.                MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
  200.      $                  ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  201.                MAXWRK = MAX( MAXWRK, M*M+M+BDSPAC )
  202.                IF( NRHS.GT.1 ) THEN
  203.                   MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
  204.                ELSE
  205.                   MAXWRK = MAX( MAXWRK, M*M+2*M )
  206.                END IF
  207.                MAXWRK = MAX( MAXWRK, M+NRHS*
  208.      $                  ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) )
  209.             ELSE
  210. *
  211. *              Path 2 - underdetermined
  212. *
  213.                MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N,
  214.      $                  -1, -1 )
  215.                MAXWRK = MAX( MAXWRK, 3*M+NRHS*
  216.      $                  ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) )
  217.                MAXWRK = MAX( MAXWRK, 3*M+M*
  218.      $                  ILAENV( 1, 'DORGBR', 'P', M, N, M, -1 ) )
  219.                MAXWRK = MAX( MAXWRK, BDSPAC )
  220.                MAXWRK = MAX( MAXWRK, N*NRHS )
  221.             END IF
  222.          END IF
  223.          MAXWRK = MAX( MINWRK, MAXWRK )
  224.          WORK( 1 ) = MAXWRK
  225.       END IF
  226. *
  227.       MINWRK = MAX( MINWRK, 1 )
  228.       IF( LWORK.LT.MINWRK )
  229.      $   INFO = -12
  230.       IF( INFO.NE.0 ) THEN
  231.          CALL XERBLA( 'DGELSS', -INFO )
  232.          RETURN
  233.       END IF
  234. *
  235. *     Quick return if possible
  236. *
  237.       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  238.          RANK = 0
  239.          RETURN
  240.       END IF
  241. *
  242. *     Get machine parameters
  243. *
  244.       EPS = DLAMCH( 'P' )
  245.       SFMIN = DLAMCH( 'S' )
  246.       SMLNUM = SFMIN / EPS
  247.       BIGNUM = ONE / SMLNUM
  248.       CALL DLABAD( SMLNUM, BIGNUM )
  249. *
  250. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  251. *
  252.       ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
  253.       IASCL = 0
  254.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  255. *
  256. *        Scale matrix norm up to SMLNUM
  257. *
  258.          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
  259.          IASCL = 1
  260.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  261. *
  262. *        Scale matrix norm down to BIGNUM
  263. *
  264.          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
  265.          IASCL = 2
  266.       ELSE IF( ANRM.EQ.ZERO ) THEN
  267. *
  268. *        Matrix all zero. Return zero solution.
  269. *
  270.          CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
  271.          CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
  272.          RANK = 0
  273.          GO TO 70
  274.       END IF
  275. *
  276. *     Scale B if max element outside range [SMLNUM,BIGNUM]
  277. *
  278.       BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
  279.       IBSCL = 0
  280.       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
  281. *
  282. *        Scale matrix norm up to SMLNUM
  283. *
  284.          CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
  285.          IBSCL = 1
  286.       ELSE IF( BNRM.GT.BIGNUM ) THEN
  287. *
  288. *        Scale matrix norm down to BIGNUM
  289. *
  290.          CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
  291.          IBSCL = 2
  292.       END IF
  293. *
  294. *     Overdetermined case
  295. *
  296.       IF( M.GE.N ) THEN
  297. *
  298. *        Path 1 - overdetermined or exactly determined
  299. *
  300.          MM = M
  301.          IF( M.GE.MNTHR ) THEN
  302. *
  303. *           Path 1a - overdetermined, with many more rows than columns
  304. *
  305.             MM = N
  306.             ITAU = 1
  307.             IWORK = ITAU + N
  308. *
  309. *           Compute A=Q*R
  310. *           (Workspace: need 2*N, prefer N+N*NB)
  311. *
  312.             CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
  313.      $                   LWORK-IWORK+1, INFO )
  314. *
  315. *           Multiply B by transpose(Q)
  316. *           (Workspace: need N+NRHS, prefer N+NRHS*NB)
  317. *
  318.             CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
  319.      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
  320. *
  321. *           Zero out below R
  322. *
  323.             IF( N.GT.1 )
  324.      $         CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
  325.          END IF
  326. *
  327.          IE = 1
  328.          ITAUQ = IE + N
  329.          ITAUP = ITAUQ + N
  330.          IWORK = ITAUP + N
  331. *
  332. *        Bidiagonalize R in A
  333. *        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
  334. *
  335.          CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
  336.      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  337.      $                INFO )
  338. *
  339. *        Multiply B by transpose of left bidiagonalizing vectors of R
  340. *        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
  341. *
  342.          CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
  343.      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
  344. *
  345. *        Generate right bidiagonalizing vectors of R in A
  346. *        (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  347. *
  348.          CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  349.      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
  350.          IWORK = IE + N
  351. *
  352. *        Perform bidiagonal QR iteration
  353. *          multiply B by transpose of left singular vectors
  354. *          compute right singular vectors in A
  355. *        (Workspace: need BDSPAC)
  356. *
  357.          CALL DBDSQR( 'U', N, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM,
  358.      $                1, B, LDB, WORK( IWORK ), INFO )
  359.          IF( INFO.NE.0 )
  360.      $      GO TO 70
  361. *
  362. *        Multiply B by reciprocals of singular values
  363. *
  364.          THR = MAX( RCOND*S( 1 ), SFMIN )
  365.          IF( RCOND.LT.ZERO )
  366.      $      THR = MAX( EPS*S( 1 ), SFMIN )
  367.          RANK = 0
  368.          DO 10 I = 1, N
  369.             IF( S( I ).GT.THR ) THEN
  370.                CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
  371.                RANK = RANK + 1
  372.             ELSE
  373.                CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
  374.             END IF
  375.    10    CONTINUE
  376. *
  377. *        Multiply B by right singular vectors
  378. *        (Workspace: need N, prefer N*NRHS)
  379. *
  380.          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
  381.             CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, A, LDA, B, LDB, ZERO,
  382.      $                  WORK, LDB )
  383.             CALL DLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
  384.          ELSE IF( NRHS.GT.1 ) THEN
  385.             CHUNK = LWORK / N
  386.             DO 20 I = 1, NRHS, CHUNK
  387.                BL = MIN( NRHS-I+1, CHUNK )
  388.                CALL DGEMM( 'T', 'N', N, BL, N, ONE, A, LDA, B, LDB,
  389.      $                     ZERO, WORK, N )
  390.                CALL DLACPY( 'G', N, BL, WORK, N, B, LDB )
  391.    20       CONTINUE
  392.          ELSE
  393.             CALL DGEMV( 'T', N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
  394.             CALL DCOPY( N, WORK, 1, B, 1 )
  395.          END IF
  396. *
  397.       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
  398.      $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
  399. *
  400. *        Path 2a - underdetermined, with many more columns than rows
  401. *        and sufficient workspace for an efficient algorithm
  402. *
  403.          LDWORK = M
  404.          IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
  405.      $       M*LDA+M+M*NRHS ) )LDWORK = LDA
  406.          ITAU = 1
  407.          IWORK = M + 1
  408. *
  409. *        Compute A=L*Q
  410. *        (Workspace: need 2*M, prefer M+M*NB)
  411. *
  412.          CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
  413.      $                LWORK-IWORK+1, INFO )
  414.          IL = IWORK
  415. *
  416. *        Copy L to WORK(IL), zeroing out above it
  417. *
  418.          CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
  419.          CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
  420.      $                LDWORK )
  421.          IE = IL + LDWORK*M
  422.          ITAUQ = IE + M
  423.          ITAUP = ITAUQ + M
  424.          IWORK = ITAUP + M
  425. *
  426. *        Bidiagonalize L in WORK(IL)
  427. *        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
  428. *
  429.          CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
  430.      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
  431.      $                LWORK-IWORK+1, INFO )
  432. *
  433. *        Multiply B by transpose of left bidiagonalizing vectors of L
  434. *        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
  435. *
  436.          CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
  437.      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
  438.      $                LWORK-IWORK+1, INFO )
  439. *
  440. *        Generate right bidiagonalizing vectors of R in WORK(IL)
  441. *        (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
  442. *
  443.          CALL DORGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
  444.      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
  445.          IWORK = IE + M
  446. *
  447. *        Perform bidiagonal QR iteration,
  448. *           computing right singular vectors of L in WORK(IL) and
  449. *           multiplying B by transpose of left singular vectors
  450. *        (Workspace: need M*M+M+BDSPAC)
  451. *
  452.          CALL DBDSQR( 'U', M, M, 0, NRHS, S, WORK( IE ), WORK( IL ),
  453.      $                LDWORK, A, LDA, B, LDB, WORK( IWORK ), INFO )
  454.          IF( INFO.NE.0 )
  455.      $      GO TO 70
  456. *
  457. *        Multiply B by reciprocals of singular values
  458. *
  459.          THR = MAX( RCOND*S( 1 ), SFMIN )
  460.          IF( RCOND.LT.ZERO )
  461.      $      THR = MAX( EPS*S( 1 ), SFMIN )
  462.          RANK = 0
  463.          DO 30 I = 1, M
  464.             IF( S( I ).GT.THR ) THEN
  465.                CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
  466.                RANK = RANK + 1
  467.             ELSE
  468.                CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
  469.             END IF
  470.    30    CONTINUE
  471.          IWORK = IE
  472. *
  473. *        Multiply B by right singular vectors of L in WORK(IL)
  474. *        (Workspace: need M*M+2*M, prefer M*M+M+M*NRHS)
  475. *
  476.          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
  477.             CALL DGEMM( 'T', 'N', M, NRHS, M, ONE, WORK( IL ), LDWORK,
  478.      $                  B, LDB, ZERO, WORK( IWORK ), LDB )
  479.             CALL DLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
  480.          ELSE IF( NRHS.GT.1 ) THEN
  481.             CHUNK = ( LWORK-IWORK+1 ) / M
  482.             DO 40 I = 1, NRHS, CHUNK
  483.                BL = MIN( NRHS-I+1, CHUNK )
  484.                CALL DGEMM( 'T', 'N', M, BL, M, ONE, WORK( IL ), LDWORK,
  485.      $                     B( 1, I ), LDB, ZERO, WORK( IWORK ), N )
  486.                CALL DLACPY( 'G', M, BL, WORK( IWORK ), N, B, LDB )
  487.    40       CONTINUE
  488.          ELSE
  489.             CALL DGEMV( 'T', M, M, ONE, WORK( IL ), LDWORK, B( 1, 1 ),
  490.      $                  1, ZERO, WORK( IWORK ), 1 )
  491.             CALL DCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
  492.          END IF
  493. *
  494. *        Zero out below first M rows of B
  495. *
  496.          CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
  497.          IWORK = ITAU + M
  498. *
  499. *        Multiply transpose(Q) by B
  500. *        (Workspace: need M+NRHS, prefer M+NRHS*NB)
  501. *
  502.          CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
  503.      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
  504. *
  505.       ELSE
  506. *
  507. *        Path 2 - remaining underdetermined cases
  508. *
  509.          IE = 1
  510.          ITAUQ = IE + M
  511.          ITAUP = ITAUQ + M
  512.          IWORK = ITAUP + M
  513. *
  514. *        Bidiagonalize A
  515. *        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
  516. *
  517.          CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
  518.      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  519.      $                INFO )
  520. *
  521. *        Multiply B by transpose of left bidiagonalizing vectors
  522. *        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
  523. *
  524.          CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
  525.      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
  526. *
  527. *        Generate right bidiagonalizing vectors in A
  528. *        (Workspace: need 4*M, prefer 3*M+M*NB)
  529. *
  530.          CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
  531.      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
  532.          IWORK = IE + M
  533. *
  534. *        Perform bidiagonal QR iteration,
  535. *           computing right singular vectors of A in A and
  536. *           multiplying B by transpose of left singular vectors
  537. *        (Workspace: need BDSPAC)
  538. *
  539.          CALL DBDSQR( 'L', M, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM,
  540.      $                1, B, LDB, WORK( IWORK ), INFO )
  541.          IF( INFO.NE.0 )
  542.      $      GO TO 70
  543. *
  544. *        Multiply B by reciprocals of singular values
  545. *
  546.          THR = MAX( RCOND*S( 1 ), SFMIN )
  547.          IF( RCOND.LT.ZERO )
  548.      $      THR = MAX( EPS*S( 1 ), SFMIN )
  549.          RANK = 0
  550.          DO 50 I = 1, M
  551.             IF( S( I ).GT.THR ) THEN
  552.                CALL DRSCL( NRHS, S( I ), B( I, 1 ), LDB )
  553.                RANK = RANK + 1
  554.             ELSE
  555.                CALL DLASET( 'F', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
  556.             END IF
  557.    50    CONTINUE
  558. *
  559. *        Multiply B by right singular vectors of A
  560. *        (Workspace: need N, prefer N*NRHS)
  561. *
  562.          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
  563.             CALL DGEMM( 'T', 'N', N, NRHS, M, ONE, A, LDA, B, LDB, ZERO,
  564.      $                  WORK, LDB )
  565.             CALL DLACPY( 'F', N, NRHS, WORK, LDB, B, LDB )
  566.          ELSE IF( NRHS.GT.1 ) THEN
  567.             CHUNK = LWORK / N
  568.             DO 60 I = 1, NRHS, CHUNK
  569.                BL = MIN( NRHS-I+1, CHUNK )
  570.                CALL DGEMM( 'T', 'N', N, BL, M, ONE, A, LDA, B( 1, I ),
  571.      $                     LDB, ZERO, WORK, N )
  572.                CALL DLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
  573.    60       CONTINUE
  574.          ELSE
  575.             CALL DGEMV( 'T', M, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
  576.             CALL DCOPY( N, WORK, 1, B, 1 )
  577.          END IF
  578.       END IF
  579. *
  580. *     Undo scaling
  581. *
  582.       IF( IASCL.EQ.1 ) THEN
  583.          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
  584.          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
  585.      $                INFO )
  586.       ELSE IF( IASCL.EQ.2 ) THEN
  587.          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
  588.          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
  589.      $                INFO )
  590.       END IF
  591.       IF( IBSCL.EQ.1 ) THEN
  592.          CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
  593.       ELSE IF( IBSCL.EQ.2 ) THEN
  594.          CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
  595.       END IF
  596. *
  597.    70 CONTINUE
  598.       WORK( 1 ) = MAXWRK
  599.       RETURN
  600. *
  601. *     End of DGELSS
  602. *
  603.       END
  604.